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1.0 


INTRODUCTION 


The ORAN program simulates a Bayesian least squares 
data reduction for orbital trajectories. It does not process 
data, but is intended to compute the accuracy of the results 
of a data reduction if measurements of a given accuracy are 
available and processed in a least squares data reduction 
program. Actual data tape input may be used, but if so, 
this provides only the time when a measurement was available 
and the estimated noise on the measurement. 

It should be noted that the ORAN program is designed 
to consider a data reduction process in which a number of 
satellite data periods are reduced simultaneously. The term 
arc refers to a specific data period over which one or more 
satellite orbits are simultaneously integrated and tracked. 

If there is more than one satellite in an arc, satellite to 
satellite tracking can be analyzed by ORAN if such tracking 
is specified on the measurement cards. 

In practice, simultaneous reduction of multiple satel- 
lite data periods would be done if parameters (such as 
station positions or geopotential coefficients) were being 
estimated which had values known to be the same for all data 
periods and all satellites analyzed. If there are no adjusted 
parameters common to all arcs, then the results for each arc 
are completely independent. 

Section 2 contains a mathematical description of the 
error analysis and Section 3 describes the partitioning of 
the error analysis equations as implemented in ORAN. Section 
4 gives a brief description of the type of measurement model 
contained in the program and Section 5 describes the force 
model equations which relate the epoch state being estimated 
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to satellite position and velocity at any 
Section 6 gives the mathematical formulae 
each force acting on the satellite and the 
variational equation used to relate errors 
parameters at epoch to errors in satellite 
velocity at any other time. 


other time, 
used to compute 
accompa nying 
in force model 
position and 
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2.0 


MATHEMATICAL DESCRIPTION 


The least squares estimator in most orbital deter- 
mination programs > assumes that N measurements can be 
modeled by the nonlinear regression equation 

— ' ff £'D * £ (2.1) 

where x contains parameters to be estimated (i.e., adjusted) 

I contains errors in parameters which are assumed to be known 
constants (i.e., unadjusted) and e is a vector of zero mean 
measurement noise. It is assumed that the covariance matrix 
associated with e is diagonal. The partitioning of parameters 
into x and x Ci.e., adjusted and unadjusted) is somewhat 
arbitrary. For any particular problem, the data will be 
insufficient to adjust all parameters subject to uncertainty 
and some reasonable subset of these parameters must be selected 
or adjustment. The final errors in the adjusted parameters 
can be decomposed into a component due to measurement noise 
and a component due to errors in the assumed values of the 
unadjusted parameters. The error statistics associated with 
the first are evaluated in the orbital determination 
program as a noise only covariance matrix. ORAN is used 

error^stati Pressing and compute 

statistics associated with the second component. 


GEODYN Program Documentation, Volume 1, Section 10. 
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Bayesian least squares estimation is characterized 
by the use of a priori information on all parameters being 
estimated. The estimation technique is otherwise identical 
to weighted least squares. In the following, no explicit 
use of a priori information will be made because of the 
complexity it introduces into the form of the equations. 
F^^I-her, when a priori information is viewed as simply 
additional observations - as is a perfectly legitimate 
procedure - then it is not necessary to treat it separately 
in the mathematical development. 

In. this section we will consider first the procedure 
used to derive the least squares estimation equations. 

The standard nonlinear least squares estimation problem is 
based on the regression equation given in (2.1). 

Since the problem must be solved iteratively using 
linear theory, for error analysis purposes the adjusted 
p 3 r ame t e r s may be ^xprossed as 


X = + 6 * ( 2 . 2 ) 

and the unadjusted parameter errors as 


- y M - Y (2.3) 


where x^ and y^ T are nominal values which are known and 
assumed to be reasonably close to the true values of x and y. 
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The linearized version of Equation (2.1) is given by 


— ^ (x^j > y^) B <$£ ~ K 6y + 


(2.4) 


where B and K are the partial derivatives of the measurement 
with respect to x and y evaluated at (x = x^, y = y^) . The 
nominal or computed measurement is defined as 


-N 


(2.5) 


by 


The weighted least squares estimate of 5x is given 


<Sx = (B T WB) -1 B T W (t_~z^j) (2.6) 


where W is the weight matrix usually assumed to be given 
by 


W" 1 = E (e e T ) 


(2.7) 


Substituting from (2.4) and (2.5) into (2.6), the 
error in estimating 6x is given by 


(6x-6x) = (B T WB) B T W(£-K6y) 


( 2 . 8 ) 


L' 
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The covariance of the estimate given by Equation (2.6) 
is defined as 

COV 6x = E[(6x-6x) (6x-<5x) T ] (2.9) 

Substituting from (2.8) into (2.9), and assuming that 
6y is uncorrelated with the measurement noise, 

E(6y e T ) = 0, 

it follows that 

COV (Sx) = (B T WB) _1 + l (B T WB) ~ 1 B T WK] COV £[(B T WB)' 1 B T WK] T (2.11) 


( 2 . 10 ) 


where COV y is the covariance matrix associated with 6y_ 


COV y = E(6y 6y T ) 


( 2 . 12 ) 


The total covariance is decomposed into a noise contribution 

and an unadjusted parameter contribution. In most practical 
♦ 

situations, the noise contribution is negligible compared to 
the unadjusted parameter effects. 

COV y is usually assumed to be a diagonal matrix 
implying that the components of 6y_ are statistically inde- 
pendent. For notational convenience in the following dis- 
cussions, the normal matrix is defined as 


N = B T W B 


(2 . 13) 


6 a/. 
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ORAN does not compute the entire systematic error 
covariance matrix given in (2.11). In order to reduce 
computational resource requirements only the diagonal 
elements of this matrix are computed. 
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SECTION 3.0 

PARTITIONING OF THE ERROR ANALYSIS EQUATIONS 


Considerable savings in both the number of computations 
and core storage requirements can be achieved by partitioning 
the matrices and vectors involved in Equations (2.6) and (2.7). 
This is done in both GEODYN and ORAN. The ORAN partitioning 
of the adjusted parameters is given by 



(3.1) 


where ou contains the ith arc orbital elements, contains 
of the other ith arc adjusted parameters, k contains 
the adjusted parameters which are common to all arcs and 
n is the total number of arcs processed. 

In a similar way the unadjusted parameter errors are 
partitioned as: 




(3.2) 
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where y^ are the unadjusted parameters associated with the 
th 1 

i arc and y^ contains the unadjusted parameters common to 
all arcs. 

It is assumed that the measurement noise is uncorrelated 
between measurements. Thus the weight matrix can be partitioned 
into arc components as : 


W 


\ 1 1 0 
# 

0 • * • w 


(3.3) 


where the individual arc weight matrices IV ^ are also diagonal. 

From the partitioning of Ax g in (3.1), the matrix of 
partial derivatives of the measurements with respect to the 
adjusted parameters can be partitioned as 


B= t B al B bl B kl 


B 0 ... 0 

B. 0 ... 0 

B, 

a l 

b l 

k l 

0 B ... 0 

OB, ... 0 

B, 

a 2 

b 2 

k 2 

0 0 . . . B 

0 0 . . . B, 

B, 

a 

b 

k 

n 

n 

n 


(3.4) 


From the partitioning of in (3.2), the matrix of partial 
derivatives of the measurements with respect to the unadjusted 
parameters can be partitioned as 
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(3.5) 


K - [KjK k ] = 


K a 0 
a l 


K 


a- 


0 0 


0 

0 


K 


n 


K 


k. 


K, 


K, 


n 


From (3.3) and (3.4) it can be shown that the normal matrix 
defined in equation (2.10) is partitioned as 


N = 


* T 

B 1 W B 

B T W B, 

B T W B. 

a a 

a b 

a 1 

B, T W B 
b a 

Bj w b, 
b b 

B b T W 

b v t w b 

k a 

b v t w b, 
k b 

B k T W B, 


(3.6) 


3.1 FIRST LEVEL PARTITIONING 

The first partitioning in ORAN separates the adjusted 
parameters into orbital elements and all other parameters as 


AX = 
— s 


a 

8 

k 


(3.1.1) 


and the corresponding partitioning of the normal matrix is 
given by 


N = 


I— 1 

n 2 

T 

N 2 

N 4 


(3.1.2) 
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where 


N, = 


N„ = 


N , = 


B T W B 


a a 


[B a T w B b 

1 B a ? W B k 

b, t w b, 
b b 

Bj w b, 
b k 

bJ w b, 
k b 

B, T W B, 
k k 


The partitioned inverse of the normal equations is given by 


N 


i = 

" M a 

M ab 


T 

M ab 

M b 


= M 


(3.1.3) 


Utilizing the relationship NM=1, where I is the identify 
matrix, it can be shown that the components of M are given 
by 


where 


M 

a 

= N x + Q M b Q 1 

(3.1.4) 

M b 

= (n 4 - n 2 t Q)' 1 

(3.1.5) 

M ab 

- *Q M b 

(3.1.6) 

Q 

"■ N l _1 N 2 



Rewriting equation (2.6) in partitioned form 

K y 


AX S 


a 



M ,t» 



T 


6 


M ab 

M b 


B 1 W 


(3.1.7) 


11 


14 < 



where 


B k - [ B k! B J 


and 


6 = 


8 

k 


Performing the matrix multiplications in (3.13) and using 
equations (3.10), (3.11) and (3.12) the errors in the orbital 
elements due to errors in unadjusted parameters are given by 


a = N i' 1 B a T W K y - Q £ 


(3.1.81 


And the errors in the rest of the adjusted parameters due to 
errors in the unadjusted parameters are given as 


B - M b [B. 1 - Q 1 B 1 ] W K y 


(3.1.9) 


3.2 SECOND LEVEL PARTITIONING 

In the second level of partitioning the common adjusted 
parameters are separated from the remaining arc adjusted 
parameters as 


8 = 



(3.2.1) 


For notational convenience define the matrix H as 


H 


I 




and note that H is a diagonal matrix 


H 



0 


0 



(3.2.2) 
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with the typical diagonal element corresponding to the ith 
arc is given by 


H i " f 1 ' \ (V T w i Ba.)' 1 B a . T W.) 


(3.2.3) 


Using the matrix partitioning shown in equations (3.6) and (3.1.2) 
equation (3.1.9) can be expanded as 



M. M, B, T W H 

b l b 2 b 


M, T M, B, T W H 

b 2 b 4 k 


^ k a —a + k n — k-^ (3.2.4) 


substituting the definitions from equation (3.8) in equation 
(3.11) it can be shown that is given by 


N. 

b l 

N, 
b 2 

-1 

M, 

b l 

M, 

b 2 


T 

N. 

b 2 

N, 

b 4 


T 

M. 

b 2 

M. 

b 4 

(3.2.5) 


where 


N, = B, 1 IV H B, 
b^ b b 


N . Bk T W H B k 
b 2 b 


I, = B, 1 W H B, 
b 4 k k 


Again using the partitioned form of the matrix inverse 


M, = N u + Q M, Q 
b-, b, x b. x 


(3.2.6) 
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(3.2.7) 



where 



Q 



N, 


2 


(3.2.8) 


Using the relationships from (3.21), (3.22) and (3.23) along with 
the partitioned forms of K and y equation (3.19) can be written 
as 

S = N bi ‘ l B b W H [K a Io * K k Y k ] - Q k (3.2.9) 

k - M b< (B k T W H - Q T B b T W H] [K a Ia * K k lk J 

3.3 ACCUMULATION OF SYSTEMATIC ERRORS AND ERROR SENSITIVITIES 

Summarizing the results of the previous two sections 
the systematic error in the common adjusted parameters is 
given by 

k = [B k T - Q T B b T ] K H [K a la * K k lk ) (3.3.1) 


14 


17 < 



The error in arc adjusted parameters (exclusive of 
orbital elements) is given by 


6 = N b ^ B b W H [K a la * K k lk ] - Q k 


(3.3.2) 


And the error in the orbital elements is given by 


a = V 1 B a T W [K a y a * K k Xk ) - Q 


(3.3.3) 


Error sensitivities are computed by taking the partial de- 
rivatives of each adjusted parameter with respect to each 
unadjusted parameter. In order to minimize core storage ORAN 
processes one arc of data at a time. Various matrices which 
involve summations over all arcs are accumulated in core. 
Components of the sensitivity matrices which are indigenous 
to individual arcs are temporarily stored on a tape or disk 
scratch file. Therefore equations (3.25), (3.26) and (3.27) 
will be used to develop these partial derivatives which will 

arc matrices . 

derivatives of the 
to the common un- 


K k (3.3.4) 


be expressed in terms of the individual 

From equation (3.25) the partial 
common adjusted parameters with respect 
adjusted parameters are given by 


3 k 


lk 


T V T T 

■ M b„ tV - <r V] w h 


In terms of the individual arc matrices can be expressed as 


M b = \t B k. w i H i r 1 - V -(V 1 w i V) 


-i 


4 *i=i 1 


i i 


-i 


(3.3.5) 


B, A W . H . ] B, 
b i 1 1 k i 
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Expanding equation (3.28) m a similar manner 





W. H. [I 



W . H . ] K 

l l J 





W. H. B 


l l 



-1 


(3.3.6) 


The partial derivatives of the common adjusted parameter with 
respect to the unadjusted parameters of the i-th arc are 
given by 



K (By 


IV. 

1 


H; B k ) 


-1 


(3.3.7) 


From equation (3.26) the partial derivatives of the i-th arc 
adjusted parameters (exclusive of the orbital element) with 
respect to the common unadjusted parameters are given by 



where 



W. H. B 


i i 



-1 



H. 

l 




3k 

(3.3.8) 


<3i ■ ( B b- 
1 


W i H i V 1 B a.' 


W. H. By 

1 1 k t 


Also from equation (3.26) the partial derivatives of the i-th 
arc adjusted parameters (exclusive of orbital elements) with 
respect to the i-th arc unadjusted parameters are given by 




W. H. B 


l l 



-1 






(3.3.9) 
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It is necessary to compute the effect of one arc's unadjusted 
parameters on another arc's adjusted parameters. Thus from 
(3.26) the partial derivatives of ith arc adjusted parameters 
(exclusive of orbital elements) with respect to the jth arc 
unadjusted parameters are given by 

9k 

-Q. (3.3.10) 

1 9 Ia. 

1 

In a similar manner from equation (3.27) the partial derivatives 
of the i-th arc orbital elements with respect to the un- 
adjusted parameters are given by 



9a. ~ 1 t ^ 93- 

- = (B 1 W. B ) B 1 W. 7K - B. - 

~ v a . l a. a. i » a. b. „ 

3l a i i 1 ( 119 ^ 

i i 


9k 


B 


- 1 3 ia. 

l 


(3.3.11) 


9a. T 1 T » 9 6.- 

- = - (B W. B ) B IV. < B, - 

~ v a . i a . a. i ) b . . 

9 — a . 1 1 1 ( 1 3 Xa 

J 

9k 


+ B 


k i 




(3.3.12) 


9a. 

—l 


9y 


-Lk 


(V w i V* B a l T “i - B b s 




3 Ik 


9k 


- B 


k. 


i 9y 


-Lk 


(3.3.13) 
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4.0 


PROPAGATION OF ERROR SENSITIVITIES 


Various ORAN output displays require- computation of 
the effects of unadjusted parameter errors on the satellite 
position and motion at some time other than epoch. This 
is given by 


Aa t 





Y 


(4.1) , 


where a is the satellite positon and velocity at time t and 
Xq is the vector of adjusted parameters, including 'the orbital 
elements at epoch time. The first term in (4.1) is due to 
errors in the estimated epoch state and the second term is 
due to the effect of the unadjusted parameters on propagation 
of that epoch state to time t. The error in the epoch state 
can be decomposed into a random component given by 


Ax r = CB T WB) ~ 1 B T W £ 


(4.2) 


and a systematic component due to errors in the unadjusted 
parameters given by 


Ax_ s = - (B T WB) _1 B T WK I 


(4.3) 


Substituting from (4.2) and (4.3) into (4.1), this 


gives 


3a 


A a. 


3 -0 


1 (B T WB) -1 B T Wc + 


3 *t 3 * t T -1 T 

— - (BMVB) 1 B iVK 

3 i 3 r 0 


X (4.4) 

2l< 


18 ^ 



The first term in (4 41 is r- a 

, ^ ' ' 15 the random error in p j 11£i 

to measurement noise. The second . -t due 

error in , i a ” ln the systematic 

-t due t0 the unadjusted parameters.- These two 

components are independent since 


E(y e 1 ) = o 


The partial derivatives anc j 9 -t 

3x„ nd 3y~ are obtai ned through 


integration of the appropriat 
described in Section S. 


e variational equations as 


The 


given by 


sensitivity of a t-n 

7 x “t t0 the unadjusted 


parameters is 


9 -t 3 £ t 


3 Y 


3 *0 


(B IVB) B f WK 


9a 


-t 


3 Y 


( 4 . 5 ) 


of each^lem^t ^ ^ CRAN t0 COmpute the sensitivity 

n in a_ t to each unadjusted parameter. The 

matrix associated with ^ is given by 


E(Aa t Aa t T ) 


3a _ 

(B T IVB)' 1 B T W 

3 -0 



C B x IVB ) " 1 B T W 


3 it 


3 Y 


3x n 


(B T 1VB) -1 b t wk 


COV Y 



9a 


3 -0 


t ^T,_ s -1 T 


(B WB)' 1 B WK 


tribution an he I” ° f (4 ' 6) iS the n0lSe con. 

parameter contribution T " ^ 

ORAN computes only the dingon^.!.™”™ 1 ^"“ ^ ‘ ' 
matrices. For each oi 0 . _ these covariance 

noise only sigma the * ^ a ^ Pr ° gram computes the 

and the totai^^a ^ *° ^*, ers 


( 4 . 6 ) 


r> O ^ 
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SECTION 5.0 

MEASUREMENT MODELING AND MEASUREMENT ERRORS 

5.1 MEASUREMENT MODELING 

In order to simulate a least squares data reduction 
ORAN must compute partial derivatives of the measurements 
with respect to both adjusted and unadjusted parameters. 

The basic types of observation in ORAN are: 

• right ascension and declination 

• range 

• range rate 

• l and m direction cosines 

• X and Y angles 

• Azimuth and elevation 

• Altimeter height 

• Inter-Satellite Range and Range Rate 

These measurements are geometric in nature. The 
computed values for the observations are obtained by applying 
geometric relationships to the computed values for the rela- 
tive positions and velocities of the satellite and the 
observer at the desired time. 


20 


23 < 



5 • 1 • 1 Range and Range Rate 
Range : 


vector : 

(5. 1.1.1) 

r is the satellite position vector (x,y,z) in 
the geocentric Earth-fixed system, and 

r ob is the station vector in the same system. 

The magnitude of this vector, p, is the (slant) 
range , which is one of the measurements . 


Consider the station- satellite 


ob 


where 


Range rate: 


The time rate of. change of this vector p is 


(S.l.1.2) 


as the velocity of the observer in the Earth-fixed sys- 
tern is zero. Let us consider that 

A 

P = pu 

(S.l.1.3) 

where 

-A 

u is the unit vector in the direction of p. 
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Thus we have 


• • s\ /\ 

p = pu + pu 


( 5 . 1 . 1 . 4 ) 


The quantity p in the above equation is the computed value 
for the range rate and is determined by 


p = u . r 


( 5 . 1 . 1 . 5 ) 


The partial derivatives of range and range rate with 
respect to the satellite position and velocity are given 
below. All are in the geocentric, Earth-fixed system. 

(The r^ refer to the Earth-fixed components of r.) 


Range : 


3p P i 

3r i p 


Range rate: 



3 P P i 

3r i p 


( 5 . 1 . 1 . 6 ) 


( 5 . 1 . 1 . 7 ') 


( 5 . 1 . 1 . 8 ) 
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The derivatives of range and range rate with 
respect to time are presented below. All are in the 
Earth-fixed system. 

Range : 



u • r 


( 5 . 1 . 1 . 9 ) 


Range Rate: 

The range rate derivative deserves special atten- 
tion. Remembering that 



( 5 . 1 . 1 . 10 ) 


We write 


P 



( 5 . 1 . 1 . 11 ) 


Thus 


•• /s • 1 •• 

p = u . p + u • p 


( 5 . 1 . 1 . 12 ) 
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Because 


P 


— (pu) = pu + pu 
dt 


(5.1.1.13) 


we may substitute in Equation 4 above for u: 


P 3 



A A II 

U • p) + U • p 


or , as 


(5.1.1.14) 


P 



(5.1.1.15) 


we may write 



(5.1.1.16) 
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The gradient of this potential with respect to the Earth- 
fixed position coordinates of the satellite is the part of 
p due to the geopotential: 




3 a e C 20 
1 


2 r 


5 sin <J>- 1-2 



r . 
* 1 


(5.1.1.17) 

We must add to this the effect of the rotation of the 
coordinate system. (The Earth-fixed coordinate system 
rotates with respect to the true of date coordinates with 
a rate 0^, the time rate of change of the Greenwich hour 
angle . ) 


The components of p are then 


p l 
p 2 
P 3 


3U . . ... 

+ [x cos 0 + y sin 01 0 + r 0 0 

3r^ 8 8 8 2 g 


(5.1.1.10) 


3U . 

+ [-x sin 8„ + y cos 0 ] 0 

3r 2 « . 8 J g 


r, 0 (5.1.1.19) 

1 g 


3U 3U 
3r^ 9 z 


(S.l.i .291 


where x and y are the true of date satellite velocity 
components . 
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5.1.2 Altimeter Height 


The altimeter height is unique in that the satellite 
is making the observation. While this is actually a 
measurement from the satellite to the surface of the Earth, 
it is taken to be a measurement of the spheroid height and 
the time rate of change of that quantity for obvious 
reasons. Using the formula for spheroid height determined 
in Section 5.1 of the GEODYN Program Documentation Volume I, 



where 

a g is the Earth's mean equatorial radius, 

f is the Earth's flattening, and 

z is r^, the z component of the Earth-fixed 
satellite vector. 
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For error analysis purposes, the partial derivatives 
of the altimeter measurement with respect to the satellite 
position, velocity, and time are needed. These are derived 
directly from the analytical expression for H ALT< 



( 5 . 1 . 2 . 2 ) 


The time derivative of altimeter range is given by 


H . 3Halt ^ . 8H ALT • + 

ALT r l — r 2 


3H 


ALT 


3r, 


3r, 


3r. 


( 5 . 1 . 2 . 3 ) 


The altimeter measurement is actually made to the 
geoid surface instead of the spheroid surface. A detailed 
geoid is necessary, however, to model the altimeter measure- 
ments to properly exploit their full accuracy. 
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5.1.3 Right Ascension and Declination 

The topocentric right ascension a and declination 
6 are inertial coordinate system measurements as illus- 
trated in Figure 4.1. ORAN computes these angles from the 
components of the Earth-fixed station-satellite vector 
and the Greenwich hour angle 0. 



The partial derivatives of these measurements with respect 
to the Earth-fixed satellite position vector r are given by 


Right Ascension: 


3a 



' p 2 



2 

+ p 2 


T 


3a ^ Pl 
3r 2 ^"l p 2 2 


(5. 1. 3. 3) 


(5. 1.3. 4) 
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Declination: 


85 

3r, 


Pi P3 


, 2 /p 2 


, + P 
1 2 



(5. 1 . 3. 6) 


(5.1.3. 7) 


(5. 1.3. 8) 


The time derivatives are given by 


Right ascension: 


a 


U 1 ^2 ~ u 2 *1 

p a-u|) 


(5. 1.3. 9) 


Declination: 


« • 

r 3 ' 0 u 3 


(S. 1. 3. 10J 


where the unit vector u is defined as 
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5.1.4 Direction Cosines 


There are three direction cosines associated with 
the station-satellite vector in the topocentric system. 
Description of these measurements requires the N, Z, and E 
(north, zenith and east baseline unit vectors which describe 
the tropocentric system along with the u) . The direction 
cosines are computed as: 


£ = u . e 

(5. 1.4.1) 

m = u . jsj 

(5. 1.4. 2) 

n = u . z 

(5. 1.4. 3) 

The Z and m direction cosines are observation types for 
ORAN . 


The partial derivatives of the direction cosines with 
respect to the satellite position vector are given by 




9r . 

l 


1 

P 



(5. 1.4. 4) 




(5. 1.4. 5) 



(5. 1.4. 6) 
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where 


- component of E in the direction 

/N 

N i = component of N in the ^ direction 

A 

= component of Z in the direction 

The time derivatives of the i and m direction cosines are 
given by 

• A « 

. p • E-Jlp 

1 = (5.1.4. 7) 

P 

• A • 

. "p • N-mp 

m = (5. 1.4. 8) 
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5.1.5 x and y Angles 


The x and y angles, as illustrated in Figure 5.2, are 
computed in a tropocentric coordinate system as 



(5. 1.5.1) 


Y = sin - '*' (m) 

a 


(5.1. 5. 2) 


The derivatives of the x and y angles with respect to the 
satellite position vector are 


3 X„ nE.-2.Z- 

a _ i i 

Sr^ p(l-m 2 ) 


(5. 1.5. 3) 



Ni -m u i 

p/l-m z 1 


(5.1. 5.4) 


and the time derivatives are given by 


4 A A 

p • (n E- 2Z) 
P (1-m 2 ) 


(5.1. 5. 5) 


Y 


a 



(5. 1. 5. 6) 
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5.1.6 Azimuth and Elevation 


Figure 5.3 illustrates the measurement of azimuth 
and elevation. These angles are computed in the topo- 
centric coordinate system as 


A 


z 


tan 


£ 


m 


(5. 1.6.1) 


E 


£ 


sin ^ (n) 


(5. 1.6. 2) 


The partial derivatives with respect to the satellite 
position vector are given by 


3 A, mE . - JIN . 

z _ 1 1 

3r i p/f-n 2 ' 


(5. 1.6. 3) 



z r pu i 

p(l-n 2 ) 


(5. 1.6. 4) 


and the partial derivatives with respect to time are 



P * (mE-£N) 
P (1-m 2 ) 


(5. 1.6. 5) 


E 


£ 



(5. 1.6. 6) 
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5 - 1 * 7 Satellite-Satellite Range and Range Rate 


The range measurement from one satellite to another 
is computed as follows. 


Let be the inertial coordinates of the transmit- 
ting satellite and X 2 the inertial coordinates of the receiving 
satellite. Then the range (or distance) between the two 
satellites is given by 



The time rate of change of range, or just range rate, 
culated by differentiating ( 1 ) with respect to time: 


R 



R 



is cal- 


CS. 1.7. 2) 


ORAN can also simulate relay range and range rate measurement 
Relay range is simply the sum of two range measurements: the 

range from some transmitting station to a satellite plus the 
range from that same satellite to another satellite. This 
configuration is given in Figure 5.4. Thus, according to the 
notation in Figure 5.4, the relay range is defined as 


s . 


R 

Relay 



(5.1. 7.3) 


Likewise, the relay range rate is the time derivation of 
(3) , or 


« 

R 

Relay 


R + R . 
1 2 


(5.1. 7.4) 
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Since the partial derivative of a sum is equal to the 
sum of the individual partial derivatives, any partial 
derivative of R Relay or R Relay can be found by summing 
the individual partials of the two quantities in the 

sum. Specifically, if one wants the partial with respect 
to some parameter a, then 


3R . 

relay 

9a 


9R 

9a 



(S.l. 7. 5) 


3R . 
relay 

Jo~ 


9R 3R 

1 

la 3a 


(S.l. 7. 6) 
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Measurement Error Parameters 


LQ 

0 

rH 

X> 

cd 

H 
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5.2 


MEASUREMENT ERRORS 


This section discusses the individual measurement 
error sources and the mathematics used to model them in 

ORAN. Table 5.1 lists all measurement error sources and 
the associated error parameters. 

5.2.1 Bias 


Bias errors are considered as constants which must 
be added on to the computed value in order to better repre- 
sent the observed. Therefore, 


Z = 
c 


Z' 

c 


+ b 


where 


(5. 2. 1.1) 


Z c is the computed measurement corrected for 

any biases 

Z c is the computed measurement based only on 
satellite geometry 

b is the bias 

Thus 

3Z 

c 

3H~ = 1 for a11 measurement types. 


(5. 2. 1.2) 



5.2.2 Timing 


Should the time tag of the measurement be incorrect, 
then a correction to this time tag is called a correction 
to timing. Any error in this correction can be found by 
computing the partial derivative of the computed measurement 
and multiplying by the time error, or 


where At is the timing error, but 


3Z 

< 

at 


az 

c 

3r 

az 

c 

3r 


3r 

at 


where 


Thus 


is the satellite position vector 
is the satellite velocity vector. 


9Z C JC 

AZ - . r At 

Tim 3F 


(5. 2. 2. 2) 


(5. 2. 2. 3) 
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5.2.3 Transit Time 


The time tag assigned to the measurement is usually 
the time at which the station receives the signal. But the 
satellite retransmitted the signal to the receiving station 
at some earlier time. Therefore, two times are involved. 

To simplify matters somewhat, the observed measurements 
usually have been corrected so the computation process of 
the computed values can assume the satellite and station at 
the same time. An error in this transit time correction is 
similar to the timing error just discussed, but now the 

systematic error is some fractional part of the estimated 
transit time, i.e. , 


AZ 

transit time 


9Z c ^ 

- 3 - • r AT 

dr 


(5. 2. 3.1) 


where AT is the error in the transit time. If p i s the 
fractional error in transit time correction, then 


AT 



(5. 2. 3. 2) 


where 


is range 

is speed of light 
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5.2.4 STATION LOCATIONS 


In the preceding sections measurement equations have 
been developed for a relative satellite-station geometry. 
These measurements are used to determine satellite position 
and motion in an inertial coordinate system at some epoch 
time. In transforming from the relative coordinate system 
at the time of measurement to the inertial system at epoch 
we must account for both the movement of the satellite and 
the movement of the station in inertial coordinates during 
the time period between measurement and epoch. The equations 
of motion for the satellite are given in Sections 5 and 6. 

The station movement is due to the movement of the Earth 
(considered as a solid body) and to the movement of the 
Earth's crust relative to the central mass. Station co- 
ordinates are referenced to a particular epoch time (usually 
1900.0) and the movement of the station since this time 
is included in the computation of a station-satellite 
measurement. The solid body component of station motion 
is due to the Earth's rotation, nutation and precession. 

These are very well known and make negligible contributions 
to station location error. 


The effect of an error in station location on 
puted measurement can be determined by the following 


the corn- 
expression 


Z 

C STA 



Ar STA 


(5. 2.4.1) 


where 


AZ, 


"STA 


is the error in the computed measurement due to 
an error in station position 


r 


is the satellite position vector 


Ar STA 


is the error in station position (in same 
coordinate system as r) 
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This is obvious if 
measurement from a 
height were raised 
would occur if the 


one considers, for example, a range 
station to a satellite. If the station 
the same effect on the measurement 
satellite height were lowered. 


At any measurement time the total station location 
error can be expressed as 


Ar STA* Ar S + Ar E + Ar SE - 4r QL ♦ 


( 5 . 2 . 4 . 2 ) 


where the components are defined as 


Ar, 


survey error. This is the error in a station's 
location relative to the local datum. Each 
station on the same local datum will have a 
different value of Ar . 


Ar T 


station location error due to uncertainty in 
location of the local datum with respect to 
the center of mass of the Earth. All stations 

on the same local datum will have the same 
Ar E . 


The remaining components are due to uncertainties in the 

movement of the Earth's crust relative to the central mass, 
lhese are: 
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Ar SE " error due to solid Earth tidal displacements, 
ihis is relatively a local error, 

Ar 0L ~ error due to ocean loading displacement. This 
is also a relatively local error which depends 
upon the distance of the station from the 
shoreline. 

Ar p = error due to polar motion 

The equations for these last three components are developed 
in the following sections. 


5.2.5 Polar Motion 


The changes in station longitude and latitude due to 
polar motion are defined as 


AX 



X 

c 


CX p SinX m " Y p COsX m ) tan <t> m (SEC ARC) 


A <t> = d> 
r m 


=-C X n C0S A + 

c pm 


V Sinx ) 

P m J 


(SEC ARC) 


(5.2.5 

(5.2.5 


where 

subscripts c and 0 denote computed and observed values 
respectively, and 

x p> Y p are angular variations (seconds of arc) in the 
position of the Earth's axis of rotation relative to 
values 3 . z. 3 . fixed epoch. 

The change in longitude, AX, can be related to UT 1 by the 
equation 
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Errors in station position due to polar motion are 

represented by uncertainties in X and Y , and are obtained 

P P 

in the form 


3 Z 3 Z 3 AA 3Z 3A<J> 

c _ c + c _ 

3X 3AA 3X 3 A 4> 3X 

P P P 


(5. 2. 5.3) 


3Z 3Z 3AA 3Z 3 A4> 

c _ c c 

3Y 3AA 3Y 3A(j> 3Y 

P P P 


(5. 2. 5.4) 


where Z c is the computed measurement. 


5.2.6 Solid Earth Tidal Displacements 

Let the total displacement of a station on the surface 
of the Earth due to the solid Earth tide be expressed in 
spherical coordinates (r,6,A): 



(5. 2.6. 1) 


N n 


r N Tr + 0 N T0 + A N ta , 


where 0 is the co-latitude 


this coordinate system the gradient 


and A if the east longitude. I n 


operator is: 


A 3 a 1 3 A 1 9 

r — + 9 + A _ 

3r r 99 r sin 0 3 A 


( 5 . 2 . 6 . 2 ) 


Let 


A 1 3 - - 1 3 

9 - — + A 

r 99 r sin 0 3 A 


( 5 . 2 . 6 . 3 ) 


Then (Diamante and Williamson, 1972): 


N n 


— r V U + — U r 

g l 2 g T 2 


(5 . 2 . 6 . 4 ) 


t K ■ ; ^ )* ? % ; 


( 5 . 2 . 6 . 5 ) 


where h £ and % 2 are ttye Love numbers of the 
kind. 


second and third 
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Using 


GM 


,, *d 2 / , " " 2 

U T = 3 r ( 3 (R^-r) - 1 , 

l 2 2 Ri 


fS.2.6.6) 


Nr 


GM d r ' 

d 3 

R d « 


/\ A 


• 3i 2 ^ R d‘ r -* 


R d + 


3 (7 - *2) <V^ 2 



rs.2.6.7) 


GM e 

Using g = — ; and r=re, the local value of the Earth's radius, 
r ^ 



(5. 2.6.8} 


N,p acts as a measurement error by introducing a time 
variation in the position of any station on the Earth's 
surface. Typically, N T is on the order of 0.5 meters, 
while the tilt d£' is on the order of 10 ** seconds of arc 
which is equivalent to horizontal displacements N~ 0.36 
meters. Values of h£ and 2 are found in the range: 
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0.587 < h 2 < 0.610 ) 

0. 068 <0.082)* (5. 2. 6. 9) 


Again, there are two contributions (S.2.6.8) 
combined effects of the Sun and Moon: 


from the 



(5.2.6.10) 


where N^, 

m 

subscript 


and N,^ are obtained from (5. 2.6. 8) replacing 
"d" with "m" and "s", respectively. 


the 
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5- 2 - 7 Ocean Loading Displacements 


Many of the currently active satellite tracking stations 
are on or near the coast. The tidal motions in the great 
bodies of water covering 70 percent of the surface of the 
Earth are now known to produce fluctuations in the observations 
of the solid Earth tidal effects. These oceanic perturbations 
account for about 10$ of the observed gravity tide, 251 of 
the observed strain tide and SOS of the observed tilt tide at 
the Earth's surface near coastlines fDiamante and Williamson, 19721 
principal ocean tide contribution to the observed surface 
gravity, strain and tilt is mainly a result of the deflection 
of the surface of the Earth under the oceanic loading. It 
is precisely this contribution to the strain tide (and hence 

station position) that is of importance in satellite opera- 
tions. 


Unfortunately, difficulties arise in developing useful 
ocean loading models of the strain tide due to: 

t 1 ) The theory of the ocean tides themselves is 

generally def icient . Complications arise from 
the fact that the oceans do not cover the entire 
surface of the planet and the depth of the oceans 
varies considerably. 

Observat ional data are limited . The state of 
the tides in the great open oceans is virtually 
unknown. The regions inland to about 400-500 km 
from continental coast-lines are dominated by 
the influence of local and regional tides. At 
about the 400 km distance from the shoreline, 
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regional tides provide about an equal contri- 
bution to the variations in the solid Earth tide 
as do the open oceans. However, further inland 
the open oceans provide the dominant contribu- 
tion. The ocean loading effect, however, drops 
off as the distance from the shore, L, in- 
creases . 

Ad ditional complications arise even where data 
are available. Even in the case of isolated 
islands, observations have been found to 
correlate with the tides of some more distant 
region, rather -than with the immediate regional 
tide. Furthermore, the most frequently avail- 
able and reliable observational data of the 
fluctuations in the solid Earth tide have been 
made with gravimeters. 

In view of these difficulties, we have developed a 
semi-empirical error model for station position which will 
incorporate available observational data. These observational 
data are mainly in the form of variations in the vertical 
component of the surface gravity which are translated into 
a model of the vertical strain variation, acceptable for 
error analysis applications. The horizontal components of 
the strain are neglected, being of secondary importance in 
any case. Kuo, et. al. (1970) have made measurements of the 
fluctuations in the M 2 and gravity tide constituents. 
However, the M 2 measurements are generally better and pro- 
vide better agreement with numerical models. Farrell (1970) 
provides some data on M 2 , 0^ , S 2 K 2 , and P-j^ gravity measure- 
ments, but the M 2 measurements are the most prevalent. 
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fraction'orrl' ”'° aSUrementS 3re presented as 

— 0f thG theoretical solid earth M, tides: 


4g t = - f 1 - + » 2 ) 


2U, 


TM, 


r - r 


( 5 . 2 . 7 .]) 


Assuming the value of the gravimetric factor T: 


'• ‘ ( 


3 

2 


i - - k 2 * h 2 


) = 


16 


( 5 . 2 . 7 . 2 ) 


and zero phase k. The result is 
value of the M 2 solid earth tide 


normalized by the theoretical 
on a rigid earth: 


6 


R 

§r 



( 5 . 2 . 7 . 3 ) 
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The measurements provide an amplitude 


^ O e ) 

in terms of a percent of 6 with a time dependence in terms 
of a phase angle 


K(r 0 ) in degrees. 


and k are indicated as functions of the local positi 


the surface of the Earth, r . 

’ e 


on on 


Kuo et. al. (1970) and Kuo and Jachens (1970) made a 
series of long-term gravimetric measurements across the con- 
tinental United States. All of the stations used in that 
study were within 1 i" of the 40" parallel of latitude. In 
employing these results, we will refer all distances from 
the coast-lines to the 40° parallel. Given the geographical 
makeup of the United States, we can expect their results to 

apply reasonably well for any U.S. station position in the 
latitude range: 

35° North < <j> < 50° North. 


Using the stations nearest to the short 


New York 


( 4 0 ° 4 9 . 0 ' N 
( 73°58.0*W 
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and 


( 38° 54 . 3 ’ N 

Point Arena, California < 

( 123°42.4'W 

as the Atlantic coast origin and Pacific coast origins (L = 0) , 
respectively, the results of Kuo, et_. al_. for the M 2 -tide 
(shown graphically in Figure 5.5) have been reduced to tabular 
form (Table 5. 2)^. Distances are measured along the 40° 
parallel . line from the Atlantic and Pacific coasts. Addi- 
tional values of the measured gravimetric factors A6 - and phase 
k for the M 2 and 0^ tides have been given by Farrell (1970) 
for a number of isolated stations, including: 

/M 2 tide: &S = 1.249, < = 5.2° 

Bermuda \ __ 

' 0 1 tide: AS = 1.198, < = 0.3° 

1.190, < = -2.0° 

1.181, < = 6.1°. 

Kuo et. al. have also made corresponding measurements for the 
0^ tidal effects across the United States. These results, 
shown graphically in Figure 5.6, do not agree as well with 
the numerical model calculations as do their M 2 tidal measure- 
ments, however. 


, M +■ i Hi* • = 

2 —a., a. 

Honolulu s _ 

(0. tide: A 8 = 


1 . 


A6 



100 


( in percent) . 
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TABLE 5.2 


Values of A3*(L) and k(L) the 


- - - , M 9 - t ide 

United Stated Referred to the 40° 


for the Continental 

and 

values (L 0) and Point Arena as the Origin 
for Pacific Coast Values (L=0)***° 



Distance 
From Pac. 
Coas t (km) 

Distance 
From Atl. 
Coas t (km) 

A3* 

K 

Degrees 

1 

~ 0 

4242 

-3.560 

3.87 

2 

202 

4040 

-0.925 

3.21 

3 

404 

3838 

0.150 

2.54 

4 

606 

3636 

0.308 

1.94 

5 

808 

3434 

0.385 

1.51 

6 

1010 

3232 

0.423 

1.21 

7 

1212 

3030 

0.308 

0.968 

8 

1414 

2.828 

0.270 

0.848 

9 

1616 

2626 

0.193 

0.786 

i n 
JL U 

*» o -t o 
lUJLU 

2424 

0.160 

0.755 

11 

2020 

2222 

0.150 

0.725 

12 

2222 

2020 

0.150 

0.725 

13 

2424 

1818 

• 0.160 

0.725 

14 

2626 

1616 

0.231 

0.725 

15 

2828 

1414 

0.308 

0.725 

16 

3030 

1212 

0.463 

0.725 

17 

3232 

1010 

0.655 

0.755 

18 

3434 

808 

1.000 

0.847 

19 

3636 

606 

1.385 

1.09 

20 

.3838 

404 

2.159 

1.45 

21 

4040 

202 

3.159 

2.18 

22 

4242 

~ 0 

4.540 

3.45 


* 

Reduced From Kuo et . al., 1970. 
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TABLE 5.3 


Values of A"? (L) and k(L) the M 2 -tide for the Continental 
United Stated Referred to the 40° Parallel of Latitude 
and with New York City as the Origin for Atlantic Coast 
Values (L=0) and Point Arena as the Origin 
for Pacific Coast Values (L=0)*** 



Distance 
From Pac. 
Coast (km) 

Distance 
From Atl. 
Coast(km) 

% 

K 

Degrees 

1 

~ 0 

4242 

-3.560 

3.87 

2 

202 

4040 

-0.925 

3.21 

3 

404 

3838 

0.150 

2.54 

4 

606 

3636 •• 

0.308 

1:94 

5 

808 

3434 

0.385 

1.51 

6 

1010 

3232 

0.423 

1.21 

7 

1212 

30 30 

0.30 8 

0.968 

8 

1414 

2828 

0.270 

0.848 

9 

1616 

2626 

0.193 

0.786 

-» r» 

XV | 

1 O “1 O 

J. U -L U 

n a A 

w'TAi-r 

0.160 

n 7c t; 

W 1 I ^ ^ 

1 

11 1 

20 20 

2222 

0.150 

0.725 

12 

2222 

2020 

0.150 

0.725 

. 13 

2424 

1818 

• 0.160 

0.725 

14 

2626 

1616 

0.231 

0.725 

15 

2820 

1414 

0.308 

0.725 

16 

3030 

1212 

0.463 

0.725 

17 

3232 

1010 

0.655 

0.755 

18 

3434 

808 

1.000 

0.847 

19 

3636 

606 

1.385 

1.09 

20 

3838 

404 

2.159 

1.45 

21 

4040 

202 

3.159 

2,18 

22 

4242 

- 0 

4.540 

3.45 


Reduced From Kuo et . al., 1970 



To good approximation the tidal potential for the M 
lunar tidal component (principal lunar semi-diurnal tide) is 
(Diamante and Williamson, 1972): 


U 


TM 


2 


3 

4 



COS <j) cos 



cos 


2t 


(5. 2. 7. 4) 


where <J> is the geocentric latitude of the station, t is the 

hour angle of the Moon and 6 is the lunar declination. d 

m m 

is the mean distance of the Moon. Let N be a unit vector 

along the north polar axis of the Earth, then: 


U, 


3 GM /r 
m 


TM, 


4 d 


2[(N-0(N-R m ) 


" * r ] 

m J 
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my 


A ^ O A A a | 

r-i r\ • r\r c ■» \ 

i. 1 " J j / 


(5. 2. 7. 5) 


Similarly, to good approximation the tidal constituent 
(principal lunar diurnal tide is: 



3 GM 
m 

4 d 


m 



sin 2<t> sin 26 cos 
m 


t . 


(5. 2. 7. 6) 
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or 




(N* r) 2 (N*R ) 2 -(N*r) (N* R ) (R • 
nr v v nr v m 



(5. 2. 7. 7) 


The measurements of Kuo et . al . , (1970) and other 
investigators are of the form: 



(5. 2. 7. 8) 


for specific lunar tidal components of U T . Equating (4) and (8) 
the solid Earth surface displacement in tlie radial direction is 
given by 


N 


SO 

Tr 




(5. 2. 7.9) 


In terms of the contributions of the M 2 and 0^ tidal effects 
(Diamante and Williamson, 1972) : 


M SO ~ V S0 , »tS0 c c\ 

N Tr ' N Tr N Tr (0 d 


4 r 2 6 x 10' 2 
__ o 

3 GM„ 


A6(M 2 ) U tm ♦ A«(0p U TQi 


r = r 


(5.2. 7. 10) 
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or 


..SO,..,. ^ ,.S0 r\i .n . M S0 rr . IX 
^Tr (-L) ^Tr ^ 2 ’ ^ ^Tr 


(5.2 


/ R. \ ** / M \ 

N~° (M ? ; L) = 1.16xl0" 2 A 6* (M,;L) [ — | ( — \ R 

\d / \ M / e 

' m ' e ' 


( 


A A A A A A 


) 


2[(N-r)(N-R*) - (R^-r)] 4 - [ 1 - ( N • R*) “ ] [ 1 - ( N * r) j [ ,(5.2 


V 


m 


n^(0 i; l) = 1.16X10' 2 Ao(O i; L) 


R. 


;4(N-?)CN*R**)[(R* A -?) - ' (N-?)(N-R**)]J , 


(5.2 


v. r ii cr? L is the distance froni the shoreline and the local phase 
angle k(L) has been taken into account by the relations: 


( R nP = CR m ) cos K (M 2 ;L) --(R m ) sin < (M 2 ;L) , ($ 2 


nr v "m y 
x x 


0**) = cy cos K ( M 2 i L ) + cy sin * CM 2 ;L) , (s .2 

y y ^ 


(R**) = (R m ) cos k (0 1 ;L) - (R) sin k (0 1 ;L) , (5.2 


.7.11) 

.7.12) 

.7.13) 

.7.14) 

.7.15) 

.7.16) 
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(R r } = ^Vv C0S K (0 1 ;L) + (R m } sin K (0 1> L) » (5.2.7.17) 

y ' x 


( R* ) = (R**) = (R ) 

v nr v in ' v nr 
z z z 


(5.2.7.18) 


For most error analysis applications, however, the phase lag 
can be assumed equal to zero and the 0^ component (equation 
3.28) may be neglected. In any case, the 0^ component can gen- 
erally be neglected for stations within the latitude band: 


-40° < <J> < 40° . 


In ORAN the 0^ component is neglected. The values of 
k and A 6 are computed by fitting a polynomial in L to 
the curves of Figure 4.5. Then Equations 12, 14, 15 
and 18 are used to compute the change in station height. 
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5.2.8 Measurement Errors Due to Medium Distortion 


In ORAN uncertainties in range and range rate due to 
the following transmission medium effects are treated. 

• Tropospheric refraction 

• Ionospheric refraction 

• Space plasma 

These are discussed in detail in the following sections. 


5.2.9 Errors in Tropospheric Refraction 

Except at very low elevation angles, (<5°) , the primary 
effect of the troposphere is a decrease in the velocity of 
propagation. At the Earth's surface, this decrease is about 
300 parts per million, decreasing to about 1 part per million 
at a height of 30 km. Considering the Earth's atmosphere to 
be horizontally stratified, as is almost always done in data 
reduction, a good approximation to the integrated tropospheric 
effect on range measurements is 


2 . 77 (N /328.5) 

AR t = meters 

.026 + sin E 


(5. 2. 9.1) 


where 


r observed r computed + AR T’ 

N is the deviation of surface index of refraction 
s 

from unity in parts per million, and 
E is the elevation angle. 
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The most serious error in applying this correction 
to data is due to errors in the surface index of refraction 
at the tracking site. For this reason, tropospheric refrac- 
tion errors are modeled in ORAN as 

3(AR t ) (2.77/328.5) (5. 2. 9. 2) 

3N g .026 + sin E 

The systematic effect of tropospheric refraction on 
range rate errors is obtained by differentiating the range 
error with respect to time. 


8(AR t ) 


3N 


(-2.77/328.5) 

j cos E E 

(.026 + sin E) 


(5. 2. 9. 3) 


Elevation : 

For elevation observations, the partial with respect 
to refraction is 


3E _ 10 3 

W ~ 16.44+930 tan E 
s 

Azimuth is not affected by refraction. 
Direction Cosines: 

3£ = -sin A sin E 3E 

3n 2 3n 

s s 


(5. 2.9. 4) 


(5. 2. 9. 5) 


3m 

3m 

s 


-cos A sin E 
z 


3E 

3m 

s 


(5. 2.9.6) 
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X and Y angles 


9X 



sin A 

2 

(sin 2 E + sin 2 A cos 2 E) 

z 


3JE 

9n 

s 


(5. 2.9. 7) 


3Y 



cos A„ sin E 
. 2 

/l - cos 2 A cos 2 E 
2 



(5. 2. 9. 8) 


5.2.10 Errors in Ionospheric Refraction 

The effect of the ionosphere on a range measurement 
is evaluated by considering 


R OBSERVED R C0MPUTED + AR I 

The correction ARj is modeled by fitting a polynomial to the 
curves in Figure 4.6, which were taken from JPL SPS 37-41, 

Volume III, page 8. The polynomial takes the form 

ARj = C Q + C 1 SIN E + C 2 SIN 2 E + C 3 SIN 3 E (5.2.10.1) 

where E is the elevation angle and the are obtained by a 
least squares fit to selected points from the curves. Errors 
in range rate are obtained by differentiating AR^. 
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RANGE CORRECTION AR T (M) 


An estimate o£ the error in ARj is given by the above 
reference to be 10% on a day-to-day basis for a particular 
location. 



Figure 5.7 Effects of Ionosphere 
on Range Measurements DSIF S-band 
System (f = 2.5 x 10^ Hz] 
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5.2.11 Errors in Space Plasma 


Space plasma represents another type of propagation 

error. Unlike the ionosphere, which is assumed to terminate 

somewhere near 600 KM above the Earth's surface, space plasma 

continues ad- inf initum , and is reasonably represented by 
2 

a 1/r law. Therefore, no closed-form solution exists for 
its effect on measurements, and an integration process must 
be performed. Let 


R 0BSERVED R C0MPUTED + AR SP 


where ARgp is modeled by the relationship 


AR 


SP 


44.3 r 

~P~ J \ is 


(5.2.11.1) 


and 


f = frequency of wave (Hz) 

Np = proton density per cubic centimeter 
s = ray path 
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PROTONS /CM 


Figure 5.7, which was taken from the JPL SPS 37-41, 

Vol. Ill, page 6, gives the assumed proton density repre- 
sentation as a function of the distance from the sun. 



Figure 5.7 Logarithmic Plot of 
Proton Density and Flux vs Distance from the Sun 

The proton density was taken to be 


N p = A[lcg R] + B 


which was converted 


into a 4 


th 


degree 


polynomial 


N = K n + K, R + K-R 2 + K,R 3 + K.R 4 


The coefficients were obtained by a least squares fit to 
selected points from Figure 5.7. 

The integration process is numerical (Simpsons Rule) 
in which the total ray path is divided in half. This 
halving process is necessary because of the possible situa- 
tion shown in Figure 5.8 in which the proton density 
curve along the ray path would be bell shaped (Figure 5.9) 



SATELLITE 


Figure 5.8 Earth-Sun-Satellite Configuration 
which Requires Halving the Ray Path 
to Accurately Represent Proton Density 
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RAY PATH (One Way) 


Figure 5.9. Proton Density Corresponding 
to Figure 5.8 


Figure 5.9 shows that a quadrature expression which does 
not consider an intermediate value of the proton density on 
the ray path would give a totally erroneous density repre- 
sentation. 
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SECTION 6.0 

FORCE MODEL AND VARIATIONAL EQUATIONS 


A fundamental part of the ORAN Program requires 
computing positions and velocities of the spacecraft at each 
observation time. The dynamics of the situation are expressed 
by the equations of motion, which provide a relationship be- 
tween the orbital elements at any given instant and the 
initial conditions of epoch. There is an additional re- 
quirement for variational partials, which are the partial 
derivatives of the instantaneous orbital elements with 
respect to the parameters at epoch. These partials are 
generated using the variational equations, which are 
analogous to the equations of motion. 


6.1 EQUATIONS OF MOTION 

In a geocentric inertial rectangular coordinate 
system, the equations of motion for a spacecraft are of 
the form 



A 


( 6 . 1 . 1 ) 


where 

r is the position vector of the satellite. 

y is the GM, where G is the gravitational constant 

and M is the mass of the Earth. 

A is the acceleration caused by the asphericity of 
the Earth, extraterrestrial gravitational forces, 
atmospheric drag, and solar radiation. 
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This provides a system of second order equations 
which, given the epoch position and velocity components, 
may be integrated to obtain the position and velocity 
at any other time. This direct integration of these 
accelerations in Cartesian coordinates is known as Cowell's 
method and is the technique used in ORAN’s orbit generator. 
This method was selected for its simplicity and its capacity 
for easily incorporating additional perturbative forces. 

There is an alternative way of expressing the above 
equations of motion: 


r = VU + . Sp + A r (6.1.2) 

where 

U is the potential field due to gravity, 

A^ contains the accelerations due to drag, and 

A^ contains the accelerations due to solar 

radiation pressure. 

This is, of course, just a regrouping of terms coupled 
with a recognition of the existence of a potential field. 

This is the form used in ORAN. 

The inertial coordinate system in which these equations 
of motion are integrated in ORAN is that system corresponding 
to the true of date system of the epoch time. The complete 
definitions for these coordinate systems (and the Earth- 
fixed system) are presented in the GEODYN Program Documen- 
tation, Volume 1, Section 3.0. 
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The evaluation of the accelerations for F is per- 
formed in the true of date system. Thus there is a require- 
ment that the inertial position and velocity output from the 
integrator be transformed to the true of date system for the 
evaluation of the accelerations, and a requirement to trans- 
form the computed accelerations from the true of date system 
to the inertial system. 


6.2 THE VARIATIONAL EQUATIONS 


The variational equations have the same relationship 

to the variational partials as the satellite position vector 

does to the equation^ of motion. The variational partials 

3x ( t 'j 

are defined as the — - — - — .where x(t) spans the true of date 

3F(t 0 ) 

position and velocity of the satellite at a given time; i.e., 


x(t) = x,y,z,x,y,z 


and SUq) spans 


x o’>V z o 


x o’>V z o 


J D 


'R 


r s 

nm’ s nm 


the epoch parameters; i.e., 

the satellite position vector at epoch 

the satellite velocity vector at epoch 

the satellite drag factor 

the satellite emissivity factor 

gravitational harmonic coefficients for 
each n, m pair 

surface density coefficients 
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GM Earth gravitational constant 

Ag lumped effect of gravitational model dif- 

ferences. The differences between several 
built-in gravity models may be propagated 
as a single parameter which approximates 
errors in low degree and order models. 
Logically, this parameter should never 
be considered adjusted. 

Yp ratio of Earth reflected radiation to 

incoming solar flux 

GM s Sun gravitational constant 

GM Moon gravitational constant 

m 

K 2 Love Number 


Also, the ORAN program can treat the perturbations in 
the Earth gravitational field due to mascon of mass m at an 
arbitrary location in the Earth. 

Let us first realize that the variational partials 
may be partitioned according to the satellite position and 
velocity vectors at the given time. Thus the required 
partials are 

aF(t) 3r(t) 

36‘(t o ) ’ 93(t 0 ) 
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where 


r(t) is the satellite position vector (x,y,z) 
in the true of date system, and 

• • • • 

r(t) is the satellite velocity vector (x,y,z) 

in the same system. 

The first of these, ■ can be obtained by the double 

30 (t ) 

integration of 0 


ar(t) 


( 6 . 2 . 1 ) 


or rathei , since the order of differentiation may be 
exchanged. 



3r(t) 


38(t 0 ) 


( 6 . 2 . 2 ) 


Note that the second set of partials, ^ may be obtained 

- 30 (t 0 ) 

by a first order integration of — . Hence we recognize 

3BCto) 3rm 

that the quantity to be integrated is _ . Using the first 

36(t 0 ) 

form given for the equations of motion in the previous 
subsection, the variational equations are given by 
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3r (t) 


9 


(6.2.3) 


9e(t o ) 3 3 (t 0 ) 


AU(t) + A R (t) + ft D (t) 


where 


U(t) is the potential field due to gravitational 
effects at time t. 

X R (t) is the acceleration due to radiation pressure 
at time t. 

A^(t) is the acceleration due to drag at time t. 


The similarity to the equations of motion is now obvious. 
When the app priate partial derivatives are computed, 
Equation (6.2.3) is of the form 


y(t) = A (t) y (t) + B (t) y(t) + C(t) 


where 


A (t) 


B (t) 


?(t) 
9r (t) 

3r(t) 
9r (t) 


(6.2.4) 
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3r(t) 

C (t) = — 

3B(t 0 ) 


3r(t) 

y(t) = — 

36(t 0 ) 

3r(t) 

y(t) = — 

33 (t Q ) 


This is a linear differential equation with coefficients A(t) , 
B(t) and C(t) which are known functions of time. In ORAN a 
Cowell predictor-corrector integrator is used to compute r(t) , 
r(t) and the coefficients A(t) , B(t) and C(t). Then an Adams- 
Cowell corrector only integrator uses these coefficients to 
solve the variational equations (6.2.4) for y(t) and y (t) . 
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SECTION 7.0 

FORCE MODEL COMPONENTS 


Equation ( 6 . 1 . 1 ) expresses the instantaneous satellite 
acceleration as the sum of individual acceleration due to 

pre e ss g ure Vlta ^ 0nal atm ° SpheriC solar radiation 

* section describes how each of these accelera- 

tors and the associated variational equations corresponding 
to Equation (6.2.4) are evaluated. Table 6.1 lists the con- 
utions to each of these acceleration components and 
the associated error parameters. 
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7.1 


THE EARTH’S POTENTIAL 


ORAN the Earth r 

spherical harmonic - Ascribed by , 

t0 m ° del loca l irregulari ties ^ lth SUrface den =ities us 




7.1.1 Spherical Harmonic Expansion 

The Earth’s potential is most conveniently expressed 
in a spherical coordinate system as is shown in Figure 7.1. 

By inspection:- 

• <J> ' , the geocentric latitude, is the angle 
measured from OX, the projection of OF in 
the X-Y plane, to the vector OF. 

• X, the east longitude, is the angle measured 
from the positive direction of the X axis 

to OS’. 

• r is the magnitude of the vector OF. 

Let us consider the point P to be the satellite 
position. Thus, OF is the geocentric Earth-fixed satellite 
vector corresponding to r, the true of date satellite 
vector, whose components are (x,y,z) . The relationship 
between the spherical coordinates (Earth-fixed) and the 
satellite position coordinates (true of date) is then 
given by 

(7.1.1) 

(7.1.2) 

(7.1.3) 
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where 9 g is the rotation angle between the true of date 
system and the Earth-fixed system (see GEODYN Program Docu- 
mentation, Volume 1, Section 3.4). 

The Earth’s gravity field is represented by the 
normal potential of an ellipsoid of revolution and small 
irregular variations, expressed by a sum of spherical 
harmonics. This formulation, used in ORAN is 



G is the universal gravitational constant, 

M is the mass of the Earth, 

r is the geocentric satellite distance, 

nmax is the upper limit for the summation (highest 
degree) , 

a g is the Earth's mean equatorial radius, 

<f> ' is the satellite geocentric latitude, 

X is the satellite east longitude, 

p^(sin<f>) indicate the associated Legendre functions, 
and 

C nm and S nm are the denormalized gravitational coefficients. 
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The gravitational accelerations in true of date co- 
ordinates (x,y,z) are computed from the geopotential, 

U(r, 4 >’,\), by the chain rule; e.g., 

3U 3r 3U 3U 3 X 

x = — — + + (7.1.5) 

3r 3x 3 4 > t 3x 3X 3x 

The accelerations y and z are determined likewise. The 
partial derivatives of U with respect to r, <p\ and X are 
given by 


(C cos mX 
' nm 


(7.1.6) 


+ S nm sin mX) (n + 1) P™ (sin 


n nmax n n 

3 U GM - — . / a A 


7 L ( 7 ) E ( S nm cos mX * C nm sin mX ) (7,1.7) 


m P n (sin 40 


3U GM nmax / a \ n n 

— , = — E h X ( bm cos raA * S nm sin mX) (7.1.6) 

n= Z m=u 


P m+ ^ (sin <t r ) - m tan <J>' P m (sin 4?) 
n n 
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The partial derivatives of r, r, and \ with respect to 
the true of date satellite position components are 





3r . 

i 


)/* 2 + y 2 ' 


dz 



3A 1 

^7" ^7 


3y y 3x 

x 3r . 
i l 


( 7 . 1 . 9 ) 


( 7 . 1 . 10 ) 


( 7 . 1 . 11 ) 
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The Legendre functions are computed via recursion 
formulae : 


Zonals : m=0 


P° (sin 40 


1 

n 


(2n-l) sin 4* 1 (sin 40 


(7.1.12) 


(n-1) P°. 2 (sin 40 


P° (sin 40 = sin <j>* 


(7.1.13) 


Tesserals: m^O and m<n 


P™ (sin 4>0 = Pjj_ 2 (sin 4»0 + (2n-l) cos <J>’ P™"* (sin 40 


P} (sin 40 


cos 4>* 


( 7 .1 .14) 

(7.1.15) 


Sectorals : m=n 


P™ = (2n- 1) cos 4>' PjJ’J (sin 40 


(7.1.16) 
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The derivative relationship is given by 


d<{> ' 


pj) (sin <f)') 


P^ +1 (sin 4>') 


m tan <£ ' P™ (sin <p ' ) 

(7.] 


It should also be noted that multiple angle formulas 
are used for evaluating the sine and cosine of mA. 


The variational equations require the computation of 
the matrix U 2 C > whose elements are given by 


U 


\ 9 2 U 

2 /iJ 9rj 


(7.] 


where 

r^ = {x, y, z}, the true of date satellite position. 
U is the geopotential. 

Because the Earth's field is in terms of r, sin <j> ' , 
and A, we write 


U 


2c 


= C 


T 9U 

sir c 


2k 


(7 


k=l 


where 


e, ranges over the elements r, sin <p ' , and A 


h 

U 7 is the matrix whose i, j element is given by 
A £ 

9 L U 


9e^ 9e. 


.17) 


. 18) 


1 . 19 ) 
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t* "h 

^ is the matrix whose i, j element is given 

■ 3 e ■ 

by 1 

3r. 

1 


and 


^2k a set three matrices whose i, j 

elements are given by 3“ e, 

Mr 


We compute the second partial derivatives of the 
potential U with respect to r, <J>', and X: 


3 2 U 2GM GM 11111 ax 


,n n 


-2 3 3 /j 

3r r r 

n=2 


— > ! Cn+1) (n+2) ^2 


m=o 


( C nm C0S mX + S nm sin mX ) P n ( sin ^ 


C 


3 2 U GM ^ x /a \ n 

= * ~2 / ( n+ l) r~) / (C cos mX 

3 r 3 4>* r 2 “ \r/ “ nm 

n= 2 v m=o 


S nm Sin mX) — (P m Csin ^ 
d <p 


(7 


3 2 u GM nmax 


n n 


3r3X r ‘ 


E c.*x, M) D 


m 


n= 2 


m=o 


( - C sin mX + S cos mX) P m (sin (Ji 1 ) (i 

nm nm n 


.1.20) 


1.21) 


1 . 22 ) 

9l- 
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where 


3 Z U GM ^ x /a \ n JL 


2 

i 

34)' 


rgB X 


(C nm cos "> A * S- 


nm 


sin mX ) 


3 

3$‘ 


( P n (= in «) 


(7.1.23) 


3 Z U GM /a \ n n 

- - E 0 £■ <■«.. ... .. 


3 4' 3 X r 


n=2 


m=o 


( 7 . 1.241 


+ 2 cos 
nm 


" A) ( P n «) 


2 

3^U 

GM 

nmax 

3V 2 " 

r 

n=2 


r 


m (C nm cos mX 


m=o 


(7.1.25) 


+ S nm sin mA ) p n ( sin <0 


3 o 


, P n (sin « 


^ nm + l 


P n (sin- <p0 - m tan <J>' p” (sin $') 


(7.1.26) 


92 < 


79 



F n m r ♦ ,|n 1 

^ \ P " tsl " +) j 


,m+ 2 


= P'“ T ^(sin <})') - (m+1) tan 0' P™ +1 (sin 00 


- m tan 0' 


P n +1 ( - sin ^ 


m tan 0* P m (sin 00 
n 


m sec 2 0 r P™ (sin 00 


(7.1.27) 


The elements of have almost been computed. 
What remains is to transform from (r, 0 1 , X) to 
(r, sin 0',X) . This affects only the partials involving 
0 ': 


3U 

3U 

30' 

3 sin 0' 

3 0' 

3 sin 0' 


(7.1.28) 


2 

3 U 


3 sin 0 


30 ’ 

3 sin 0* 



30' 


3U 


2 

3 0 ' 


3 sin 0 1 3 0' 3 sin 0 1 


7Z 


where 


(7.1.29) 


3 0' 

= sec 0' 

3 sin 0' 


(7.1.30) 


3 2 0 ' 


3 sin 0' 


2 


sin 0' sec -5 0' 


(7.1.31) 
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For the and C ^ matrices, the partials of r, 
sin <t>', and X are obtained from the usual formulas: 



We have for C^: 

3r r. 

x 


3 sin <J>* ' -zr. 1 dz 



(7.1.32) 


(7.1.33) 


(7.1.34) 


(7.1.35) 


(7.1.36) 


(7.1.37) 
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The C 2k are s y mme ' tI 'ic . The necessary elements 
are given by 


3 2 r 


3 r. 3 r . 
i J 


ILL i . i !li 


r 3 r. 

J 


(7.1.38) 


3 2 sin <J>' 

3r. 3r. 
i J 


3 2 X 


3r. 3r. 
i J 


3 z r. r. 
i L 


- 2r. 


(x 2 +y 2 ) 2 


2 2 
x +y 


~T 

r 


ay 


3 r. 

l 


r. 

J 


3 z 


3 . 

l 


3 z 

+ r. + z 

l 


3x 


3 r. 
l 


3x 3y 


ay 3 x 


3r. 3r. 

J J 


dr. dr. 
J J 


3 . 

J 


dr. 

l 


J 


(7.1.39) 


(7.1.40) 


If gravitational constants, C nm or S nm are being 
treated, we require their partials in the f matrix 

for the variational equations computations. These 
partials are 


3C 


nra 



cos (nu) P^Csin (7 j 41) 


3C 


nm 


GM /a 


n 


- — I- m — I — ] sin (mX) P (sin q 1 ) 
ax/ r \ r / n 


(7.1.42) 
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ac 


nm 


3U\ GM /aV 1 

_ , - - M ‘O* (»*) 


P n +1 ( sin ^ 


m tan <p' P™ (sin 4>') 


(7.1.43) 


The partials for S nm are identical with cos (mX) re- 
placed by sin (mX) and with sin (mX) replaced by 
-cos (mX ) . 

These partials are converted to inertial true of 
date coordinates using the chain rule; e.g.. 


3C 


nm 


/ 3U\ _ _ 

( 3xJ 3C 


3u\ 3r 


u i o r 3 /- 3 u\ 3 X 

7 3X 3C nm \ 3 7 Sx 


r -3U\ 34>' 


3C nm \ 3< *>7 3 X 


(7.1.44) 
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The logic and appropriate array dimensions of the 
ORAN program have been modified to accommodate uncertainties 
in the coefficients of the zonal and tesseral harmonic terms 
in the geopotential up to and inclduing (30,30). This is 
near to the practical limit for the current 360 series 
machines. Beyond (30,30), the accuracy of the computations 
of the perturbations and partials becomes poor for double 
precision and excessive amounts of computing time are involved 
in calculating the Legendre polynomials. However, the 
highly localized gravitational anomalies that require higher 
order spherical harmoncis for adequate representation can 
be simulated by using the surface density model described 
in Section 7.2. 
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7.2 


MASCONS AND SURFACE LAYER DENSITY 


The surface layer density represents the acceleration 
due to a localized gravity anomaly. This acceleration is pro- 
portional to the area of the particular block in question and 
inversely proportional to the square of the distance from the 
spacecraft. The surface layer density block can be repre- 
sented as a mascon by equating the errors in the spacecraft 
acceleration due to a mascon and local density layers. Let 

Rimtx- be this acceleration. Then: 

DEN 


r den. 

3 



(7.2.1) 


where <j>. is the density of block i; A. is the area of block i; 

^ 1 3 

x^ is the distance along the j axis; and is the cube of 

distance from block i to the spacecraft. 

ORAN was modified to account for both propagations 
of errors in this model and adiustments of the local density 

V * 

values. Any number of blocks can be generated by the program. 
However, the user is cautioned not to specify so many blocks 
that an excessive amount of core is required. 

For integration of these force model partials, the 
differential equation 



SF 

3x 

3F 


— 

34> i 

+ 

(7.2.2) 

3x 

H j. 


must be solved. 


To do this, ORAN required the second term on the 
right-hand side. It is found by differentiating with 

respect to 4^ , 
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A. X. . 
i n 


d<p . 


The complete details of the surface layer density model in the 
ORAN program have been given by Martin (1972) . 

MASCONS are modeled as point masses which may be located 
at any point in the earth. These point masses perturb the satel- 
lite orbit in the same manner as do geopotential harmonics and 
are consequently treated similarly in the ORAN program. The 
gravitational potential of a MASCON produces a satellite force 
which can be integrated numerically along with the variational 
equations for epoch element and geopotential coefficient partial 
derivatives. Thus, as far as the EASST program is concerned, 
the only thing unique about a MASCON is its forcing function in 
the variational equations. The appropriate forcing function will 
now be derived. 

The potential of a mass m at a point which is a distance 
p from the mass is, by definition 


where G is the gravitational constant, 
of a mass (M-m) at a distance is 


Similarly, the potential 


G (M-m) 
P i 


The potential of both masses at the same point, P, is 


y _ G (M-m ) + Gm 


(7.2.6) 
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This expression for the potential V can be used to 
derive the disturbing potential of a MASCON. Let m be the 
mass of the MASCON and (M-m) be the mass of the body containing 
the MASCON, so that the total mass within the primary body is 
M. Then the potential of the central force is GM/r where r 
is the distance of the point P from the center of gravity of 
the MASCON and the primary body. If the MASCON (m) is fixed 
with respect to the mass (M-m) , then we can write the disturbing 
potential as 

AV = V - GM/r 


_ G (M-m) + Gm GM 

" p i p " r 


(7.2.7) 


This expression can be simplified if we are prepared 
to neglect terms of orderf^j . Even for the moonl^j is of 
order 10 s . For the earth, for which/^,) is much smaller, we 
are completely justified in neglecting J^j . With this ap- 
proximation the disturbing potential can be written as: 


AV 


Gm 

P 


— (xX + yY + zZ) 
r 3 


Gm 

r 


(7.2.8) 


where (x,y,z) are the coordinates of P, and (X,Y,Z) are the 
coordinates of the MASCON relative to the center of gravity of 
the system. 


From the disturbing potential we can derive the compo- 
nents of the disturbing acceleration by differentiation. The 
x-component is 


3 (AV) 
3x 


Gm 

„ 3 


(x - X) 


Gm ( 
r 3 \ 


(x-X) 


— (xX + yY + 
r 2 



(7.2.9) 


with identical forms for the other two components. 
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7.3 


LUNAR AND SOLAR PERTURBATIONS 


Consider a system of n point masses m. (i = 0 , 1 , 2 , . . . n- 1) . 
The equations of motion of the i mass relative to the mass 
uIq can be written: 


Oi 


G(m 0 +m.) 

~ — r„ . + G 

„ •> Oi 

r 0i 


n-3 




j = l 
j r i 


r . . 

il 




r J r ^ 

ij 0j 


(7.3.1) 


Let m Q denote the Earth and nu the satellite. The satellite's 
motion relative to the Earth- is then (m. is negligible with 
respect to the other masses) 





+ 



r ij _ r o; 


(7.3.2) 


The variational equation for an error in the gravitational 
constant Gm^ (j = 2,3,...n-l) is then 


3r 


01 


5r 01 
3 (Gnu ) 


3 CGm. ) y 


dt‘ 


111 

r l3 


11 

‘0j 


(7.3.3) 


These variational equations are integrated numerically to 
give 3(r 01 )/3(Gmj) . 
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Recent determinations of GM for the Moon give a value of 
about 4902.6 kg-Vsec^ with a standard deviation of about 0.25 
kg'Vsec^. Consequently a reasonable estimate of the uncertainty 
in GM for the Moon is given by 0.25/4902.6 ~ 5.0x10'^. 
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7.4 


SOLID EARTH TIDES 


. „ The tldal displacement of mass produces a perturbation 

1,1 th<S field of the Earth, U F ol lo „i„ g the 

small response assumption. Love (13n) represented the dis- 
turbance potential, at the surface r=R , by 


w * k mv 

(7.4.1) 


This can be generalized (Kaula, 1968) as 


w ■ E 


n=2 


k n W 


(7.4.2) 


Since K e << 1 
approximation, 


terms involving n>2 can be neglected to first 
and 


W •* k 2 W’ 


(7.4.3) 


There 

value 

based 


is an unresolved discrepancy with the best theoretical 
of k 2 which is equal to 0.290 and the best estimates 
on satellite measurements, which now center on 0.25. 
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The second order third body disturbing potential at 
the surface of the earth is given by 


GM, R 
d e 


(cos 0) 


where G is the gravitational constant, is the mass of the 
disturbing body, R g is the radius of the Earth, and R^ is 
the magnitude of the vector R^ from the Earth to the disturbing 
body. The angle 9 is the angle enclosed by R d and the vector 
r from the earth to the satellite. Evaluating the Legrendre 
polynomial and using the vector formulation for cos 0 yields 


GM, R 
d e 


— 1 2 
r 


R d r 


Letting the 


notation refer to the unit vector. 


k ? GM, R * ^ 7 

U n (r) = — -4 - 4 - 3 (R d *r) Z - 1 

u 2 R d r |_ 


The acceleration produced at position r by the solid 
Earth tidal bulge is (Diamante and Williamson, 1972): 


r = VU r 


3 GM, R" 
1 d e 

7 k 2 73 4 

2 R d r 




A A \ /s 

3 (R d * r ) **- 1 J r 


(R d -r) 
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The variational equation is: 


A A A 


d 2 /3r \ 1 GMj R 5 . * 

77 VlkTy = 7 7“ 7 ’ [3 - 15( V r) lr * 6 (r h'^, 

N L* a 


k, GM, R 
2 a e 


/\ A 


U 3 r 


1 0 5 ( R d • r ) r- 


3k, 


A A / A d r 

• 3 °( R d* r M R d- 


3k 


2 


A 8r M 

-A 

15 

r + 

V »v. 




3 r 


3k. 


r A 9 T 

30 (R , *r) [r 

1 ' 3k, 


R 


A A A 

37 ) 

3 - 15 (R, »r ) i 

V 

U. 

3 k 2 j 



V F 


(7.4.8) 


The variational equation is actually given by the sum 
of the expression above for the Moon and Sun, i.e., with the 
subscript "d" replaced by "s" and "m". These equations were 
derived under the assumption that the response time of the Earth 
is negligible. Actually, the tide will be displaced by some 
angle 6 = n Q • At where n is the angular velocity of the Earth 
and At is the time lag of the tide. This phase lag may be 
accounted for by replacing R^ in the above Equation with a vector R^ 

which corresponds to rotated through an angle 6 in the direc- 
tion of Earth's rotation (Diamante and Williamson, 1972). 
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7.5 


LUMPED GEOPOTENTIAL ERROR DUE TO GRAVITY MODEL 
DIFFERENCES 


The ORAN program possesses the capability of modeling 
the total set of errors in the set of spherical harmonic 
coefficients used to perform an orbital data reduction. 

The rationale for the use of this form is basically as 
follows. Let represent a satellite orbit generated 

or estimated using a set of geopotential coefficients 
estimated by an investigator, or group of investigators, 
which we have denoted by APL. Similarly, let x g^Q represent 
the same orbit ephemeris obtained (estimated) in the same 
way except that geopotential coefficients estimated by a 
group denoted by SAO were used. Let us further assume that 
group APL and group SAO operated completely independently 
using different tracking data on different satellites. 

On this basis, we deduce that errors in their geopotential 
models and thus geopotential dependent errors in x^p^ and 
x SAO are ^dependent . It is not necessary to assume that 
the two models are equally accurate. 


Let us then consider the variance of the difference 

We can obtain this variance in terms 
of the variance of the individual model errors by first 
writing the difference as 


between x A p L and x SA q. 


X APL ‘ X SAO ( X APL " X P ' fx SAO“ x t^ 


6x APL ' 6x SAO (7.5.1) 


where 


is the true ephemeris and Sx^p^ and <$Xg A Q are the 
errors in the two orbits using the two different geo- 
potential models. Then 
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(7.5.2) 


V ar ^ X APL -X SA(P ^ t ( x apt _x c;An^ ^ x APr~ x QAn) 3 


APL SACr *■ APL A SAO' 


Var 6x ApL + Var <5x SAQ 


For the special case in which the models are of equal variance 
(accuracy) , we have 


Var (<5x ApL ) = Var (<5x SA0 ) = j Var ( x 


APL~ X SAO J 


(7.5.3) 


If the models are not of equal accuracy, but we can relate 
the variances of the two models by a constant factor, 


Var 6x ApL k Var <5x SAQ , 
we then have the relation 


Var («x SA0 ) = jij Var Cx ApL -x SA0 ) 


If the SAO model is the more accurate, then 


(7.5.4) 


(7.5.5) 


1 

rnr 


< 


i 

i 


(7.5.6) 


and we have an upper limit on the coefficient of the gravity 
model difference variance which is needed to obtain the vari- 
ance of the SAO model. 
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In practice, we would not expect k to be exactly a 
constant, nor would we expect to find geopotential models 
with completely independent errors, but we would expect to 
be able to choose models which are nearly independent and 
with a suitable definition of the variances of an d 

SXgAo C e * g - » sampled over all possible orbits), then k is 
a constant. In any event, we will use the above relation 
obtained using the SAO gravity model, in terms of the gravity 
model difference effect. 

It is possible to compute the effects on an orbit 
of a gravity model difference using the integration of a 
set of variational equations. Consider first the variational 
equation for a spherical harmonic which we can write as 



• • •• 

3r(t) 3r(t) + 3r(t) 

3r(t) 3C nra 3r(t) 


3r(t) + 3r(t) 


3C 


nm 


3C_ 


nm 


( 7 . 5 . 6 ) 


Where r(t) is the satellite position (x,y,z) at time t and 
r(t) is the acceleration. Note that the acceleration due to 
gravity is dependent only on position. Therefore, the second 
term on the right hand side of the above equation is zero. 

The effect of a small error 6C on the orbital parameters 
at time t is given by 




!l + 3 1 I 

ac ac ( 

nm nm / 


6C 

nm 


( 7 . 5 . 7 ) 


A corresponding equation holds for the S coefficient of 

nm 

the same degree and order 
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(Sr(t)) = 


(7.5.8) 



I 9 r(t) 3r(t) + 3r(t) ) 6S 

) 3r(t) as as ( nm 

v v 7 nm nm ) 


If we sum (7.5.7) an ^ ( 7 , 5 , 8 ) for a range of values of n and m 
we get the variational equation for the total error in the 

orbital parameters at time t due to errors in the spherical 
harmonic coefficients 



(6r(t) ) 


E 

n, m 


3r(t) 

3r(t) 


3r(t) &r 3r(t) 



+ 


3r(t) 

^nm 


5S 


nm 


+ 


3r(t) 

3C 


nm 


6C 


nm 


(7.5.9) 


Within the limits of linearity (and ORAN is a linear error 
analysis program) , the quantities in brackets are the orbit 
differences due to a set of differences (or errors) in the 
spherical harmonic coefficients. We specialize this set of 
differences to be obtained by differencing geopotential 
model sets which are both as accurate and independent as 
possible. Equation (7.5.9) then gives the variational equation 
for the effects which the set of coefficient differences pro- 
duce on the estimated orbit. 


In ORAN each spherical harmonic coefficient may be 
1) adjusted, 2) treated as an individual unadjusted parameter 
or 3) included in the lumped parameter model of equation (7.5.9). 
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7.6 


DIRECT SOLAR RADIATION PRESSURE 


For a spherical satellite with an integrated reflectivity 
pL , the acceleration is (Diamante, 1972b): 


_ - - S ( k \ 2 ( R S- r ’ 

a = r = -v(l + p R ) - - d t 

R c\mj ° |R s -F|" 


where 

0, when the satellite is in the Earth's 
shadow , 

1, when the satellite is illuminated by 
Sun, 

r = satellite position vector (Earth centered 
coordinate system), 

Eg = position vector of the sun (Earth centered 
coordinate system) , 

A = cross-sectional area of the satellite^, 
m = mass of the satellite. 



The cross-sectional area is in a plane normal to the direc- 
tion of the radiation flux. A is therefore constant for 
spherical satellites. However, for other geometries, vari- 
able cross-section along the orbit must be taken into 
account or an average value of A adopted. 
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Clearly, when >_ 0, the satellite is always in full sun- 
light. If f^ < 0, the satellite will be in shadow when: 


r *r 


— ^ 2 2 
(r*R s )' £ 


(7.6. 2) 


R is the radius of the Earth, Let: 
e 


? * — 2 ? 
*2 ’ r ' <V r) - R e ' 


Then 


1 , £ 1 > 0 


v = \ 1, < 0 and f ^ > 0 


\C , f^ < 0 and f£ £ 0 


(7.6.3) 


(7.6.4) 


An error in the direct solar radiation pressure force 
model can be expressed as an uncertainty in p^, or as an 
uncertainty in C^. Then 





3r 


1 + C 


R 


3C 


R 


V 


C R S- r > 
I R S' r I 


+ 



(7.6.5) 
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where V is the gradient operator and F includes the effects 
of all the other forces acting on the satellite. The vari- 

3 x* 

ational equation for -jp— reduces to (Diamante, 1972b): 


3 » d‘ 


3C 


R 


dt‘ 


3r 


,3C 


R 


= -v 


-U 

,m/ 



1 + 3 


C R (Fs-r) 3r 


R s -r| 2 3 C r 


'R 


|R s -r| 3 


3r 

3C 


R 


3r 


3C 


F . 


R 


(7.6.C) 


7.7 EARTH REFLECTED RADIATION 

Over the long term, the radiation budget of the Earth 
must be in equilibrium. A large fraction of the solar radia- 
tion received by the Earth is reflected back, almost immedi- 
ately, at nearly the same wavelengths at which it is received 
The ratio of the reflected radiation to the incoming solar 
flux is known as the albedo , . The best estimate of the 

long-term mean global albedo for the Earth has been obtained 
from TIROS VII measurements: 


R 


0.32 


(7.7.1) 



The remaining fraction of the incoming solar flux is absorbed 
by the atmosphere, oceans and surface of the Earth and is 
eventually re-emitted as the far-infrared planetary radiation. 
From equilibrium requirements, this fraction is: 


y IR = 1 " y R = 0,68 (7.7.2) 

Part of the reflected radiation is scattered diffusely 
and part is scattered specularly. The specular reflection is 
mostly from small bodies of water scattered over the surface 
of the Earth and occasional calm areas of the oceans . There- 
fore , 


Y 


R 


D 


+ 


Y S ’ 


(7.7.3) 


where the subscripts D and S stand for diffuse and spectral. 
Using a value of 4/3 for the index of refraction of water, 
the refractivity at normal incidence is 0.02. 

Furthermore, although 71% of the Earth's surface is covered 
by water, only a small fraction of that can be expected to 
produce specular reflection at any one instant in time. 
Therefore, to good approximation, 


Y d * Y r = 0.32 (7.7.4) 

Since for the purposes of error analysis we need only 
consider the radial accelerations, we have employed a simplified 
analytic model for the Earth reflected radiation (Diamante, 

1972b) . 
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Danjon (1954) made observations of the Earth-shine 
reflected from the dark hemisphere of the Moon. He found 
that the phase function of the Earth is better represented 

by a non-Lambert Law expression. At great distances, the 

• 2 
Danjon relation can be approximated by a (l-cos8 s ) dependence. 

The Earth receives solar energy at a rate: 


ttR 


e 


2 



ergs/ sec. 


and reflects: 



ergs/sec 


away. If this is distributed at great distances according 
*to £ unc £ i. o n i 


C (r ) (l-cos8 s ) 2 , 


then the total energy flux through a spherical surface of 
radius r, assuming that the flux is radially outward, is: 


IT 

2irr 2 C(r) J* sin s (l-2 cosB s + cos 2 8 s )d6 s 
0 


16 


ir r 


C(r) . 
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Setting this equal 


to 


y d"V 



( 7 . 7 . 5 ) 


Therefore, the Wyatt-Danjon expression for the radiation 
pressure is 



( 7 . 7 . 6 ) 


and the radial acceleration is: 



\ S 


3/16 Y d ( 1 + P r ) - 


m c 



(1-cosB r. 
s 


( 7 . 7 . 7 ) 


A force error in the Earth-reflected radiation pressure 
can be modeled as an error in Yq* Noting that: 


A A 

cosB s = R g • r 


( 7 . 7 . 8 ) 


where the notation again denotes the unit vector, the 

variational equation is then: 
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or 


3r 




(a r ^ 


3 

— C1+P R ) 
16 R 


A S 
m c 



+ 




(7.7.10) 


where 



(7.7.11) 
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7.8 


ATMOSPHERIC DRAG 


The acceleration r on a satellite due to the atmos- 
pheric drag force is given by the expression: 


-1/2 C 


D 



r 


r 


(7.8.1) 


where 

Cp is the coefficient of drag for satellite, 

A Q is the cross-section area of the satellite, 

(m ) 

mg is the mass of the satellite (kg) 

r is the velocity vector of the satellite, 
relative to the atmosphere (under the 
assumption that the atmosphere rotates 
with the Earth) (m/sec) , 

r is the position vector from the center of 
the Earth to the satellite (m) , 

p is the mass density of the atmosphere 

3 

at position r (kg/m ) . 


1 ? 8 


104 



Strictly speaking, Ag is the area of the satellite projected 

in the V direction. However, in practice the mean cross- 

• r 

section is used (constant) and Cp is also taken to be a 
satellite dependent constant. 

The major source of error in determining the drag 
accelerations arises from uncertainties in the density, p. 
Since these uncertainties can usually be expected to exceed 
drag coefficient uncertainties, density errors can con- 
veniently be represented as errors in C^. Let 


F 


2 \m s / 


r 


r 


(7.8.2) 


The variational equations for 



are 


then : 


3r 


3C 


D 


7? 



F + 




F + 




F + 


3r 




F 



(7.8.3) 


where F involves the effects of all other forces (assumed 
independent of velocity) and V is the gradient operator. 
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The righthand side of equation (2.3) may be expanded, 
remembering that p is not a function of velocity. The ex- 
pressions for the partial derivatives have been given by 
Williamson, Martin and Dutcher, 1971. All derivatives of 
p, however, must be evaluated numerically for a given satel- 
lite orbit using the atmospheric density model described 
below. An error in the satellite position due to an uncer- 
tainty in Cp (drag error) is then obtained by integrating 

3 IT 

the. above variational equations for . 

The Jacchia 1965 Model Atmosphere, with revisions 
through 1968, has been selected to represent the atmospheric 
density in the drag error calculations of the ORAN program. 

The atmospheric density calculation in the ORAN pro- 
gram is performed by subroutine DENSTY. This subroutine is 

a modification of the one currently employed in the GEODYN 
programs and is based on the static diffusion models of 

Jacchia (1964, 1968). 

The 1968 Jacchia Model Atmosphere allows the calcula- 
tion to be structured into two major steps. In the first 
step, the exospheric temperature is calculated from data on 
solar flux, geomagnetic index, solar ephemeris and satellite 
position information. In the second step, the density is 
calculated as a function of the exospheric temperature T , 
and the satellite altitude above the surface of the Earth, z. 
This second step is accomplished by using a polynomial fit to 
tabulated solutions of differential equations. The complete 
procedure is summarized in Figure 2.1. 
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The daily value of the 10.7 cm solar flux strength, 

F 10 y, is in units of 10 watts/m /cycle/second of band- 
width. In the GEODYN program, the daily values of F^q y 
are data input and F^g y is obtained by averaging the daily 
values over the preceding 54 days (two solar rotations) . 

In the ORAN program, it is anticipated that for most error 
analysis applications and mission planning, the daily values 
will either be unavailable (before the fact) or represent an 
unnecessary burden on storage and computation for the level 
of accuracy required. Consequently, in the ORAN program, 

% 

F 10 . 7 = F 10.7- 

Average values of F^g y have been obtained by averaging the 
daily record of F^g y over intervals of the 11-year solar 
cycle corresponding to periods of high, medium and low solar 
activity. These values are 218.80, 130.62 and 124.35, re- 
spectively. Commonly accepted values of F^q y for very high, 
medium and very low levels of solar activity are 250-275, 

150, and 65-75, respectively. Appropriate values may be 
selected for the anticipated solar activity levels in the 
period of interest. 

The correction for geomagnetic activity was provided 
by Jacchia, Slowey and Verniani, 1967 and is based on ap , 
the three-hour planetary geomagnetic index. In applying this 
result in the ORAN program, the three-hour equivalent plane- 
tary amplitude ap was replaced by the daily equivalent plane- 
tary AP in order to minimize data input and storage. The 
index Ap is defined as the arithmetic average of the eight 
values of ap for a given day. For application of the 
geomagnetic correction in ORAN, the same remarks apply here 
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as for input data on F^q y . Average values of Ap corres- 
ponding to the same averaging periods over the 11 year solar 
cycle employed in the calculation of average values of y 

for periods of high, medium and low solar activity were 
found. These are respectively Ap=20.24, 13.63 and 10.98, 
corresponding to the values F 1Q y=218.80, 130.62 and 124.35 
for high, medium and low solar activity. The indices ap and 
Ap can range from 0 to 400. Days when Ap<2 are considered 
to be geomagnetical ly quiet and slightly disturbed days are 
those for 4<Ap<^10 . Moderate disturbances are in the ran^e 
Ap=40-50 and intense disturbances range Ap_100. 
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7.9 


MICROMETEORITE IMPACT AND THRUST ERRORS 


Collision of a spacecraft with micrometeorites can 
be modeled by uncertainties in spacecraft thrusting. If 
is the acceleration due to the thrusting, then the error 
in position and velocity due to an uncertainty in r T at the 
time of application is 


Ar 


3r 



1 

- (Ar T ) 

2 1 


(At) 2 , 


» 3 x •• 

Ar = Ar T = (ir^,) (At) 

3 x* fp 


where At is the duration of the thrust and Ar T is the uncerta 
in the magnitude. These relations provide the initial values 
for integrating variational equations to provide the effects 
the orbit as a function of time. 

Information on the meteoric environment near to the 
Earth and in interplanetary space, useful for application of 
the ORAN error model, is provided in Space and Planetary 
Environment Criteria Guidelines For Use In Space Vehicle 
Development, 1971 Revision. 


(7.9.1) 

(7.9.2) 

inty 
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